﻿// 2024/03/12 待实现以下第2段Strassen代码。
// 链接： https://blog.csdn.net/weixin_52042488/article/details/126345027 


namespace ConsoleAlgorithms.Utils
{
    public class StrassenMatrixMultiplication: ITestable
    {
        public void Test()
        {
            Random random = new Random();
            int i = 0;
            int[] a = new int[size * size];
            int[] b = new int[size * size];
            int[] c = new int[size * size];
            //时间随机数种子生成  
            //给A和B矩阵赋随机值 （0~10） 
            for (i = 0; i < (size * size); i++)
            {
                a[i] = random.Next(10);
                b[i] = random.Next(10);
            }
            //打印矩阵A和B 
            Console.WriteLine("\nA matrix is:");
            for (i = 0; i < (size * size); i++)
            {
                Console.Write($"{a[i],5}");
                if (((i + 1) % size) == 0)
                    Console.WriteLine(); 
            }

            Console.WriteLine("\nB matrix is:");
            for (i = 0; i < (size * size); i++)
            {
                Console.Write($"{b[i],5}");
                if (((i + 1) % size) == 0)
                    Console.WriteLine();
            }

            //进行矩阵乘法 
            Matrix_Multiply(a, b, c, 0, 0, size);
            //打印结果矩阵C 
            Console.WriteLine("\nC matrix is:");
            for (i = 0; i < (size * size); i++)
            {
                Console.Write($"{c[i],5}");
                if (((i + 1) % size) == 0)
                    Console.WriteLine();
            }
            Console.WriteLine();
        }

        //宏定义改变矩阵的大小
        //此处size指的是n*n矩阵的宽度或者长度n
        //由于Strassen算法本身的局限性
        //两个相乘的矩阵只能是n*n,且n只能是2的幂
        //即size为1,2,4,8,16...
        int size = 4;


        //矩阵的合并
        //就是将A11,A12,A21,A22合并为
        //	|A11	A12| 
        // 	|A21	A22|
        //这样的形式 
        void Merge_Matrix(int[] a, int[] b, int[] c, int[] d, int[] c0, int rows)
        {
            //rows是子矩阵的宽度，那么合并后矩阵的宽度就是2*rows 
            int i = 0;//子矩阵遍历索引 
            int j = 0;//合并矩阵遍历索引，此处先合并A11 
            for (i = 0; i < (rows * rows); i++)
            {
                //如果执行了rows次就需要换行，要锁定到合并后矩阵的第二行，所以加上2*rows即可 
                if ((i % rows == 0) && i != 0)
                {
                    j = (rows * 2) * (i / rows);//(i/rows)代表需要换多少行；(rows*2)就是行数 
                    c0[j] = a[i];
                    j++;
                }
                else
                {
                    c0[j] = a[i];
                    j++;
                }
            }
            //此处是A12的首元素，令索引等于子矩阵的那个函数
            j = rows;
            for (i = 0; i < (rows * rows); i++)
            {
                //这里的条件和之前的一样的，其它四部分也是一样的，底层逻辑是一致的 
                if ((i % rows == 0) && i != 0)
                {
                    j = (rows * 2) * (i / rows) + rows;//换行要换到第二部分，所以需要再加一个rows 
                    c0[j] = b[i];
                    j++;
                }
                else
                {
                    c0[j] = b[i];
                    j++;
                }
            }
            //此处是A21的首元素，令索引等于子矩阵的那个函数 
            j = rows * 2 * rows;
            for (i = 0; i < (rows * rows); i++)
            {
                if ((i % rows == 0) && i != 0)
                {
                    j = rows * 2 * rows + (rows * 2) * (i / rows);
                    c0[j] = c[i];
                    j++;
                }
                else
                {
                    c0[j] = c[i];
                    j++;
                }
            }
            //此处是A22的首元素，令索引等于子矩阵的那个函数
            j = rows * 2 * rows + rows;
            for (i = 0; i < (rows * rows); i++)
            {
                if ((i % rows == 0) && i != 0)
                {
                    j = rows * 2 * rows + rows + (rows * 2) * (i / rows);
                    c0[j] = d[i];
                    j++;
                }
                else
                {
                    c0[j] = d[i];
                    j++;
                }
            }
        }


        //两个大小相同的矩阵的加法 
        void Matrix_Add(int[] x, int[] y, int[] c0, int rows)
        {
            int i = 0;
            for (i = 0; i < (rows * rows); i++)
            {
                c0[i] = x[i] + y[i];
            }
        }


        //矩阵的乘法递归函数 
        void Matrix_Multiply(int[] x, int[] y, int[] c0, int loca, int locb, int rows)
        {
            //进行索引的申明
            //递归过程为了节约时间，没有进行原矩阵的切割赋值
            //直接采用索引的方法进行矩阵的伪切割 
            int loca11, loca12, loca21, loca22;
            int locb11, locb12, locb21, locb22;
            //结果矩阵的四分子矩阵申明 
            int[] c11 = new int[(rows / 2) * (rows / 2)];
            int[] c12 = new int[(rows / 2) * (rows / 2)];
            int[] c21 = new int[(rows / 2) * (rows / 2)];
            int[] c22 = new int[(rows / 2) * (rows / 2)];
            //中间暂存矩阵申明 
            int[] temp1 = new int[(rows / 2) * (rows / 2)];
            int[] temp2 = new int[(rows / 2) * (rows / 2)];

            //四分矩阵后的子矩阵宽度 
            int newrows = rows / 2;
            //如果矩阵的宽度为1，那么就直接元素相乘，递归结束 
            if (rows == 1)
            {
                c0[0] = x[(loca)] * y[(locb)];
            }
            else
            {
                //索引的计算，进行矩阵伪分割 
                //因为n为2的幂，所以直接二分即可 
                loca11 = loca;
                locb11 = locb;
                loca12 = loca + (rows / 2);
                locb12 = locb + (rows / 2);
                loca21 = loca + (size * (rows / 2));
                locb21 = locb + (size * (rows / 2));
                loca22 = loca + (size * (rows / 2)) + (rows / 2);
                locb22 = locb + (size * (rows / 2)) + (rows / 2);
                //根据算法公式进行运算
                //涉及矩阵相乘的进入递归 
                Matrix_Multiply(x, y, temp1, loca11, locb11, newrows);
                Matrix_Multiply(x, y, temp2, loca12, locb21, newrows);
                Matrix_Add(temp1, temp2, c11, newrows);
                Matrix_Multiply(x, y, temp1, loca11, locb12, newrows);
                Matrix_Multiply(x, y, temp2, loca12, locb22, newrows);
                Matrix_Add(temp1, temp2, c12, newrows);
                Matrix_Multiply(x, y, temp1, loca21, locb11, newrows);
                Matrix_Multiply(x, y, temp2, loca22, locb21, newrows);
                Matrix_Add(temp1, temp2, c21, newrows);
                Matrix_Multiply(x, y, temp1, loca21, locb12, newrows);
                Matrix_Multiply(x, y, temp2, loca22, locb22, newrows);
                Matrix_Add(temp1, temp2, c22, newrows);
                //递归后，进行矩阵的合并 
                Merge_Matrix(c11, c12, c21, c22, c0, newrows);
            }
        }

    }
}
